library(runjags)
library(MASS)
library(parallel)

###Load data
load("pmrdata.RData")

##########################
###########Table 1########

##########Proportion fluent in Russian
1-mean(pmrdata$rf,na.rm=T)
1-mean(pmrdata$rf[pmrdata$Russian==1],na.rm=T)
1-mean(pmrdata$rf[pmrdata$Moldovan==1],na.rm=T)
1-mean(pmrdata$rf[pmrdata$Ukrainian==1],na.rm=T)

##########Proportion fluent in Moldovan
mean(pmrdata$mf,na.rm=T)
mean(pmrdata$mf[pmrdata$Russian==1],na.rm=T)
mean(pmrdata$mf[pmrdata$Moldovan==1],na.rm=T)
mean(pmrdata$mf[pmrdata$Ukrainian==1],na.rm=T)

##########Proportion fluent in Ukrainian
mean(pmrdata$uf,na.rm=T)
mean(pmrdata$uf[pmrdata$Russian==1],na.rm=T)
mean(pmrdata$uf[pmrdata$Moldovan==1],na.rm=T)
mean(pmrdata$uf[pmrdata$Ukrainian==1],na.rm=T)

################Sums of different groups
nrow(pmrdata)
sum(pmrdata$Russian,na.rm=T)
sum(pmrdata$Moldovan,na.rm=T)
sum(pmrdata$Ukrainian,na.rm=T)



###########################################################
###########################################################
#Main experimental and observational analyses; analyses of#
#NAs#######################################################
###Figures 3 and 5, Appendices I, J, K##################### 
###########################################################

####Experimental balance checks
summary(lm(tg=="Control" ~ rf + mf + uf + NRussian + Ukrainian + Other + income + 
             age + male + highered + urban, data=pmrdata))
summary(lm(tg=="T1" ~ rf + mf + uf + NRussian + Ukrainian + Other + income + 
             age + male + highered + urban, data=pmrdata))
summary(lm(tg=="T2" ~ rf + mf + uf + NRussian + Ukrainian + Other + income + 
             age + male + highered + urban, data=pmrdata))
summary(lm(tg=="T3" ~ rf + mf + uf + NRussian + Ukrainian + Other + income + 
             age + male + highered + urban, data=pmrdata))
####No consistent evidence of balance problems

####Threshold initial values for experimental outcome
taustar<-polr(method="probit", as.factor(DV1) ~ rf + mf + uf + NRussian + Ukrainian + Other + 
                income + age + male + highered + urban, data=pmrdata)$zeta

###Observational threshold initial values
gammastar<-matrix(nrow=10,ncol=4)
for(i in 1:10){
  gammastar[i,]<-polr(method="probit", as.factor(pmrdata[,i+23]) ~ pmrdata$rf + pmrdata$mf + 
                        pmrdata$uf + pmrdata$NRussian + pmrdata$Ukrainian + pmrdata$Other + 
                        pmrdata$income + pmrdata$age + pmrdata$male + pmrdata$highered + 
                        pmrdata$urban)$zeta
}
gammastar<-apply(gammastar,2,mean)

###########Data file
hmod<-list(N=nrow(pmrdata), tg=pmrdata$tg,
           y1=pmrdata$DV1, y2=as.matrix(pmrdata[,24:33]), 
           y1NA=pmrdata$DV1na, y2NA=as.matrix(pmrdata[,34:43]), 
           g0=rep(0,4), G0=diag(.01,4), 
           t0=rep(0,3), T0=diag(.01,3),
           rf=pmrdata$rf, mf=pmrdata$mf, uf=pmrdata$uf, NRussian=pmrdata$NRussian,
           Other=pmrdata$Other, Ukrainian=pmrdata$Ukrainian, Male=pmrdata$male, 
           Age=pmrdata$age, Urban=pmrdata$urban,
           NIOTA=12, I0=diag(1,12), i0=rep(0,12), 
           NALPHA=7, A0=diag(1,7), a0=rep(0,7), 
           B0=diag(1,5), b0=rep(0,5),
           HigherEd=pmrdata$highered, Married=pmrdata$married, Locality=pmrdata$locality,
           Prof=pmrdata$prof, Resp=pmrdata$resp, House=pmrdata$house, ymin=pmrdata$ymin, 
           ymax=pmrdata$ymax,
           bI0=rep(0,5), BI0=diag(1,5), ae1=rep(1,14))

##hmodel <- run.jags(method="parallel",model="pF2.txt", inits=list(gammastar=gammastar,
##                                               taustar=taustar), 
##                   monitor=c("alpha", "beta", "tau", "mu.a", "iota", "gamma", "alphaNA", 
##                             "betaNA", "mu.aNA", "iotaNA", "pr", "pm",  "pu", "bI", 
##                             "betaI","gammaI","alphaI", "pIe", "pIm","e1", "rIl","rIh",
##                             "Income"), data=hmod, n.chains=4, adapt=5000, burnin=50000, 
##                   sample=1000, thin=500, summarise=F,plots=F, modules=c("glm", "lecuyer"))

##save(hmodel, file = "pF2.RData")
load("pF2.RData")


library(rjags)
library(ggplot2)

#####################Check convergence
h1<-as.mcmc.list(hmodel)
gd<-gelman.diag(h1[,1:395],multivariate = F) ###Exclude imputed income bc many values at ceiling or floor
sum(gd[[1]][,1] > 1.1) 
####All model parameters have rhat < 1.1

#####################Convert to single chain
cds<-as.mcmc(hmodel)

#####Estimate income and age average across posterior draws
inc<- apply(as.matrix(cds[,396:972]),1,mean)
age<-mean(pmrdata$age)

############################experimental analysis
####Function to estimate probability respondent supports PMR independence across MCMC draws 
#with different ethnic and linguistic characteristics, holding everything else constant at mean
#or mode
bpred<-function(cont,t1,t2,t3,eps,beta,inc,age){
  mu<-matrix(nrow=nrow(cont),ncol=36)
  pred<-matrix(nrow=nrow(cont),ncol=36)
  for(i in 1:nrow(cont)){
    ###monolingual Russians
    mu[i,1] <-cont[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,2] <-t1[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,3] <-t2[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,4] <-t3[i,]%*%rbind(1,0,0,0,0,0)
    ###bilingual R/M Russians
    mu[i,5] <-cont[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,6] <-t1[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,7] <-t2[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,8] <-t3[i,]%*%rbind(1,0,1,0,0,0)
    ###bilingual R/U Russians
    mu[i,9] <-cont[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,10] <-t1[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,11] <-t2[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,12] <-t3[i,]%*%rbind(1,0,0,1,0,0)
    ###monolingual R Moldovans
    mu[i,13] <-cont[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,14] <-t1[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,15] <-t2[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,16] <-t3[i,]%*%rbind(1,0,0,0,1,0)
    ###bilingual R/M Moldovans
    mu[i,17] <-cont[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,18] <-t1[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,19] <-t2[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,20] <-t3[i,]%*%rbind(1,0,1,0,1,0)
    ###monolingual M Moldovans 
    mu[i,21] <-cont[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,22] <-t1[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,23] <-t2[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,24] <-t3[i,]%*%rbind(1,1,1,0,1,0)
    ###monolingual R Ukrainians
    mu[i,25] <-cont[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,26] <-t1[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,27] <-t2[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,28] <-t3[i,]%*%rbind(1,0,0,0,1,1)
    ###bilingual R/U Ukrainians
    mu[i,29] <-cont[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,30] <-t1[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,31] <-t2[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,32] <-t3[i,]%*%rbind(1,0,0,1,1,1)
    ###monolingual U Ukrainians
    mu[i,33] <-cont[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,34] <-t1[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,35] <-t2[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,36] <-t3[i,]%*%rbind(1,1,0,1,1,1)
    for(j in 1:36){
      pred[i,j]<-  1- pnorm(eps[i]  - (mu[i,j] + beta[i,]%*%rbind(inc[i],age)))
    }
  }
  return(pred)
}

###Third threshold
eps<-as.matrix(cds[,36])
####Parameters by experimental conditions
cont<- as.matrix(cds[,seq(1,24,4)])
t1<- as.matrix(cds[,seq(2,24,4)])
t2<- as.matrix(cds[,seq(3,24,4)])
t3<- as.matrix(cds[,seq(4,24,4)])
####Treatment-invariant parameters
beta<- as.matrix(cds[,29:30])

###Apply function
mu<-bpred(cont,t1,t2,t3,eps,beta,inc,age)



##########Table H.1
###Russian prime on monolingual Russian-speaking Russian
c(median(mu[,3]-mu[,1]), HPDinterval(as.mcmc(mu[,3]-mu[,1]),prob=.95))
###Moldovan prime on bilingual Russian/Moldovan Russians, also in txt
c(median(mu[,6]-mu[,5]), HPDinterval(as.mcmc(mu[,6]-mu[,5]),prob=.95))
###Ukrainian prime on bilingual Russian/Ukrainian Russians
c(median(mu[,12]-mu[,9]), HPDinterval(as.mcmc(mu[,12]-mu[,9]),prob=.95))

###Moldovan prime on monolingual R Moldovans
c(median(mu[,14]-mu[,13]), HPDinterval(as.mcmc(mu[,14]-mu[,13]),prob=.95))
###Moldovan prime on bilingual R/M Moldovans, also in txt
c(median(mu[,18]-mu[,17]), HPDinterval(as.mcmc(mu[,18]-mu[,17]),prob=.95))
###Moldovan prime on monolingual M Moldovans, also in txt
c(median(mu[,22]-mu[,21]), HPDinterval(as.mcmc(mu[,22]-mu[,21]),prob=.95))

###Ukrainian prime on monolingual R Ukrainians
c(median(mu[,28]-mu[,25]), HPDinterval(as.mcmc(mu[,28]-mu[,25]),prob=.95))
###Ukrainian prime on bilingual R/U Ukrainians
c(median(mu[,32]-mu[,29]), HPDinterval(as.mcmc(mu[,32]-mu[,29]),prob=.95))
###Ukrainian prime on monolingual U Ukrainians
c(median(mu[,36]-mu[,33]), HPDinterval(as.mcmc(mu[,36]-mu[,33]),prob=.95))

###Estimate posterior median and hpd
mum<-apply(mu,2,median)
muhpd<-HPDinterval(as.mcmc(mu),prob=.95)

#####Create dataset 
Condition<-rep(c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian"),9)
Ethnicity<-c(rep("Russian",12),rep("Moldovan",12),rep("Ukrainian",12))
Language<-factor(c(rep("Russian monolingual",4),rep("Russian/Moldovan bilingual",4),
            rep("Russian/Ukrainian bilingual",4), rep("Russian monolingual",4),
            rep("Russian/Moldovan bilingual",4), rep("Moldovan monolingual",4),
            rep("Russian monolingual",4),rep("Russian/Ukrainian bilingual",4),
            rep("Ukrainian monolingual",4)))
mu<-data.frame(mum,muhpd,Condition, Ethnicity,Language)
names(mu)[1]<-"value"

####Control condition graphic
muc<-mu[mu$Condition=="Control",]
muc$levels<-interaction(muc$Ethnicity,muc$Language)
muc$levels<-droplevels(muc$levels)
muc$levels = factor(muc$levels,levels(muc$levels)[c(3,7,6,4,8,9,2,5,1)])
levels(muc$levels)<-c("  Russian monolingual", "  Russian/Ukrainian bilingual", 
                      "  Russian/Moldovan bilingual", " Russian monolingual",        
                      " Russian/Ukrainian bilingual", " Ukrainian monolingual",      
                      "Russian monolingual", "Russian/Moldovan bilingual",  
                      "Moldovan monolingual")
muc$Ethnicity<-as.character(muc$Ethnicity)

####Figure J.1
ggplot(data = muc, aes(levels,value,color=Ethnicity, shape=Ethnicity))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(3.5,6.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="") + 
  theme(text = element_text(size=15), legend.position=c(.25,.25))  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

###subgroup analyses
mu$levels<-interaction(mu$Ethnicity,mu$Language,mu$Condition)
mu$levels<-droplevels(mu$levels)
mu$levels = factor(mu$levels,sort(levels(mu$levels)))
mu$levels = factor(mu$levels,levels(mu$levels)[c(13:16,21:24,17:20,25:36,5:12,1:4)])
levels(mu$levels)<-c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian", 
                     " Control", " T_Moldovan", " T_Russian", " T_Ukrainian",
                     "  Control", "  T_Moldovan", "  T_Russian", "  T_Ukrainian",
                     "   Control", "   T_Moldovan", "   T_Russian", "   T_Ukrainian",
                     "    Control", "    T_Moldovan", "    T_Russian", "    T_Ukrainian",
                     "     Control", "     T_Moldovan", "     T_Russian", "     T_Ukrainian",
                     "      Control", "      T_Moldovan", "      T_Russian", "      T_Ukrainian",
                     "       Control", "       T_Moldovan", "       T_Russian", "       T_Ukrainian",
                     "        Control", "        T_Moldovan", "        T_Russian", "        T_Ukrainian")
mu$Ethnicity<-as.character(mu$Ethnicity)

#########Figure 4. Russian
ggplot(data = mu[which(mu$Ethnicity=="Russian" & mu$Condition != "T_Ukrainian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(1.5,3.5,6.5)) + 
  geom_segment(aes(y=mu[9,1], yend=mu[9,1], x=3.5, xend=6.5)) +
  geom_segment(aes(y=mu[1,1], yend=mu[1,1], x=.5, xend=3.5)) +
  geom_segment(aes(y=mu[5,1], yend=mu[5,1], x=6.5, xend=9.5)) +
  theme_bw() + labs(y="", x="", color="",shape="") + guides(colour = guide_legend(nrow = 2))+ 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

#############Figure 4, Moldovan
ggplot(data = mu[which(mu$Ethnicity=="Moldovan" & mu$Condition != "T_Ukrainian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(1.5,3.5,6.5)) + 
  geom_segment(aes(y=mu[17,1], yend=mu[17,1], x=3.5, xend=6.5)) +
  geom_segment(aes(y=mu[13,1], yend=mu[13,1], x=.5, xend=3.5)) +
  geom_segment(aes(y=mu[21,1], yend=mu[21,1], x=6.5, xend=9.5)) +
  theme_bw() + labs(y="", x="", color="",shape="") + guides(colour = guide_legend(nrow = 2))+ 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

##########Figure J.2
ggplot(data = mu[which(mu$Ethnicity=="Ukrainian"),], 
       aes(levels,value,shape=Language,  color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(1.5,4.5,5.5,8.5,9.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="", linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

#########Figure J.4
ggplot(data = mu[which(mu$Ethnicity=="Russian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(1.5,4.5,5.5,8.5,9.5)) + 
  geom_segment(aes(y=mu[9,1], yend=mu[9,1], x=4.5, xend=8.5)) +
  geom_segment(aes(y=mu[1,1], yend=mu[1,1], x=.5, xend=4.5)) +
  geom_segment(aes(y=mu[5,1], yend=mu[5,1], x=8.5, xend=12.5)) +
  theme_bw() + labs(y="", x="", color="",shape="") + guides(colour = guide_legend(nrow = 2))+ 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

#############Figure J.3
ggplot(data = mu[which(mu$Ethnicity=="Moldovan"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(1.5,4.5,5.5,8.5,9.5)) + 
  geom_segment(aes(y=mu[17,1], yend=mu[17,1], x=4.5, xend=8.5)) +
  geom_segment(aes(y=mu[13,1], yend=mu[13,1], x=.5, xend=4.5)) +
  geom_segment(aes(y=mu[21,1], yend=mu[21,1], x=8.5, xend=12.5)) +
  theme_bw() + labs(y="", x="", color="",shape="") + guides(colour = guide_legend(nrow = 2))+ 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

###################################
############Table J.3##############
#Manually edit to correctly format#
###################################

library(stargazer)
pm<-apply(cds,2,median)
hpd<-HPDinterval(cds, prob=.9)

##########Experimental coefficients
pes<-cbind(pm[seq(1,28,4)],1, hpd[seq(1,28,4),1], 2,hpd[seq(1,28,4),2], 3,
           pm[seq(2,28,4)],1, hpd[seq(2,28,4),1], 2,hpd[seq(2,28,4),2], 3,
           pm[seq(3,28,4)],1, hpd[seq(3,28,4),1], 2,hpd[seq(3,28,4),2], 3,
           pm[seq(4,28,4)],1, hpd[seq(4,28,4),1], 2,hpd[seq(4,28,4),2], 4)
rownames(pes)<-c("Intercept", "Not fluent in Russian", "Fluent in Moldovan", 
                 "Fluent in Ukrainian", "Not Russian", "Ukrainian", 
                 "Other")
stargazer(pes, digits=2)

################Condition-invariant coefficients
pes<-cbind(pm[29:36],1, hpd[29:36,1], 2,hpd[29:36,2], 4)
rownames(pes)<-c("ln(Income)", "ln(Age)", "Male", 
                 "Higher education", "Urban", "gamma_1",
                 "gamma_2","gamma_3")
stargazer(pes, digits=2)

###############################################################################
##################################Observational analyses#######################
###############################################################################

####Function to estimate probability respondent supports given separatist outcomes across MCMC 
#draws with different ethnic and linguistic characteristics, holding everything else constant at 
#mean or mode
bpreds<-function(cont,eps,inc,age){
  mu<-matrix(nrow=nrow(cont),ncol=9)
  pred<-matrix(nrow=nrow(cont),ncol=9)
  for(i in 1:nrow(cont)){
    ###monolingual Russians
    mu[i,1] <-cont[i,]%*%rbind(1,0,0,0,0,0,inc[i],age)
    ###bilingual R/M Russians
    mu[i,2] <-cont[i,]%*%rbind(1,0,1,0,0,0,inc[i],age)
    ###bilingual R/U Russians
    mu[i,3] <-cont[i,]%*%rbind(1,0,0,1,0,0,inc[i],age)
    ###monolingual R Moldovans
    mu[i,4] <-cont[i,]%*%rbind(1,0,0,0,1,0,inc[i],age)
    ###bilingual R/M Moldovans
    mu[i,5] <-cont[i,]%*%rbind(1,0,1,0,1,0,inc[i],age)
    ###monolingual M Moldovans 
    mu[i,6] <-cont[i,]%*%rbind(1,1,1,0,1,0,inc[i],age)
    ###monolingual R Ukrainians
    mu[i,7] <-cont[i,]%*%rbind(1,0,0,0,1,1,inc[i],age)
    ###bilingual R/U Ukrainians
    mu[i,8] <-cont[i,]%*%rbind(1,0,0,1,1,1,inc[i],age)
    ###monolingual U Ukrainians
    mu[i,9] <-cont[i,]%*%rbind(1,1,0,1,1,1,inc[i],age)
    for(j in 1:9){
      pred[i,j]<-  1- pnorm(eps[i]  - mu[i,j])
    }
  }
  return(pred)
}

####Third threshold (neither/nor vs support)
eps<-as.matrix(cds[,166])
####################Integration w Moldova, Figure I.1
contMol<- as.matrix(cds[,c(seq(44,103,10),seq(114,133,10))])
####################Autonomy w/in Moldova, Figure I.1
contAut<- as.matrix(cds[,c(seq(45,103,10),seq(115,133,10))])
####################Confederation w Moldova, Figure I.1
contCon<- as.matrix(cds[,c(seq(46,103,10),seq(116,133,10))])
####################Independence, Figure 2
contInd<- as.matrix(cds[,c(seq(47,103,10),seq(117,133,10))])
####################Integration w Russia, Figure I.2
contRus<- as.matrix(cds[,c(seq(48,103,10),seq(118,133,10))])
####################Autonomy w/in Russia, Figure I.2
contRA<- as.matrix(cds[,c(seq(49,103,10),seq(119,133,10))])
####################Confederation w Russia, Figure I.2
contRC<- as.matrix(cds[,c(seq(50,103,10),seq(120,133,10))])
####################Integration w Ukraine, Figure I.2
contUkr<- as.matrix(cds[,c(seq(51,103,10),seq(121,133,10))])
####################Autonomy w/in Ukraine, Figure I.3
contUA<- as.matrix(cds[,c(seq(52,103,10),seq(122,133,10))])
####################Confederation w Ukraine, Figure I.3
contUC<- as.matrix(cds[,c(seq(53,103,10),seq(123,133,10))])

##########################Figure 2, Table H.1, txt stats####################################
#To create figures/ estimate stats for different outcomes, replace contInd w relevant matrix
############################################################################################
mu<-bpreds(contInd,eps,inc,age)

##################Table H.1, txt
#############Difference btw bilingual R/M Russian and monolingual R Russian
c(median(mu[,2]-mu[,1]), HPDinterval(as.mcmc(mu[,2]-mu[,1]),prob=.95))
#############Difference btw bilingual R/U Russian and monolingual R Russian
c(median(mu[,3]-mu[,1]), HPDinterval(as.mcmc(mu[,3]-mu[,1]),prob=.95))
#############Difference btw monolingual R Moldovan and Russian, in txt
c(median(mu[,4]-mu[,1]), HPDinterval(as.mcmc(mu[,4]-mu[,1]),prob=.95))
#############Difference btw bilingual Moldovan and monolingual R Moldovan
c(median(mu[,5]-mu[,4]), HPDinterval(as.mcmc(mu[,5]-mu[,4]),prob=.95))
#############Difference btw monolingual M Moldovan and monolingual R Moldovan, in txt
c(median(mu[,6]-mu[,4]), HPDinterval(as.mcmc(mu[,6]-mu[,4]),prob=.95))
#############Difference btw monolingual R Ukrainian and Russian
c(median(mu[,7]-mu[,1]), HPDinterval(as.mcmc(mu[,7]-mu[,1]),prob=.95))
#############Difference btw bilingual R/U and monolingual R Ukrainian
c(median(mu[,8]-mu[,7]), HPDinterval(as.mcmc(mu[,8]-mu[,7]),prob=.95))
#############Difference btw monolingual U and monolingual R Ukrainian, in txt
c(median(mu[,9]-mu[,7]), HPDinterval(as.mcmc(mu[,9]-mu[,7]),prob=.95))

mum<-apply(mu,2,median)
muhpd<-HPDinterval(as.mcmc(mu),prob=.95)


Language<-as.factor(c("Russian monolingual","Russian/Moldovan bilingual","Russian/Ukrainian bilingual",
                      "Russian monolingual", "Russian/Moldovan bilingual", "Moldovan monolingual",
                      "Russian monolingual", "Russian/Ukrainian bilingual", "Ukrainian monolingual"))

Ethnicity <- as.factor(c(rep("Russian", 3),rep("Moldovan", 3),rep("Ukrainian", 3)))

mu<-data.frame(mum,muhpd,Language,Ethnicity)

names(mu)[1]<-"value"

mu$levels<-interaction(mu$Ethnicity,mu$Language)
mu$levels<-droplevels(mu$levels)
mu$levels = factor(mu$levels,levels(mu$levels)[c(3,7,6,4,8,9,2,5,1)])
levels(mu$levels)<-c("  Russian monolingual", "  Russian/Ukrainian bilingual", 
                     "  Russian/Moldovan bilingual", " Russian monolingual",        
                     " Russian/Ukrainian bilingual", " Ukrainian monolingual",      
                     "Russian monolingual", "Russian/Moldovan bilingual",  
                     "Moldovan monolingual")

mu$Ethnicity<-as.character(mu$Ethnicity)


ggplot(data = mu, aes(levels,value,color=Ethnicity, shape=Ethnicity))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(3.5,6.5)) + 
  theme_bw() + labs(x="", y="", color="",shape="") + 
  theme(text = element_text(size=15), legend.position=c(.25,.25))  + 
  coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17)) 



################################Table I.2#######################
###########Edit manually to create tables#######################
################################################################

################Integration w Moldova
pes<-cbind(pm[seq(44,163,10)], 1, hpd[seq(44,163,10),1], 2, hpd[seq(44,163,10),2], 3,
           pm[seq(45,163,10)], 1, hpd[seq(45,163,10),1], 2, hpd[seq(45,163,10),2], 3,
           pm[seq(46,163,10)], 1, hpd[seq(46,163,10),1], 2, hpd[seq(46,163,10),2], 3,
           pm[seq(47,163,10)], 1, hpd[seq(47,163,10),1], 2, hpd[seq(47,163,10),2], 4)
rownames(pes)<-c("Intercept", "Not fluent in Russian", "Fluent in Moldovan", 
                 "Fluent in Ukrainian", "Not Russian", "Ukrainian", 
                 "Other","ln(Income)", "ln(Age)", "Male", 
                 "Higher education", "Urban")
stargazer(pes, digits=2)

################Integration w Russia
pes<-cbind(pm[seq(48,163,10)], 1, hpd[seq(48,163,10),1], 2, hpd[seq(48,163,10),2], 3,
           pm[seq(49,163,10)], 1, hpd[seq(49,163,10),1], 2, hpd[seq(49,163,10),2], 3,
           pm[seq(50,163,10)], 1, hpd[seq(50,163,10),1], 2, hpd[seq(50,163,10),2], 4)
rownames(pes)<-c("Intercept", "Not fluent in Russian", "Fluent in Moldovan", 
                 "Fluent in Ukrainian", "Not Russian", "Ukrainian", 
                 "Other","ln(Income)", "ln(Age)", "Male", 
                 "Higher education", "Urban")
stargazer(pes, digits=2)

################Integration w Ukraine
pes<-cbind(pm[seq(51,163,10)], 1, hpd[seq(51,163,10),1], 2, hpd[seq(51,163,10),2], 3,
           pm[seq(52,163,10)], 1, hpd[seq(52,163,10),1], 2, hpd[seq(52,163,10),2], 3,
           pm[seq(53,163,10)], 1, hpd[seq(53,163,10),1], 2, hpd[seq(53,163,10),2], 4)
rownames(pes)<-c("Intercept", "Not fluent in Russian", "Fluent in Moldovan", 
                 "Fluent in Ukrainian", "Not Russian", "Ukrainian", 
                 "Other","ln(Income)", "ln(Age)", "Male", 
                 "Higher education", "Urban")
stargazer(pes, digits=2)

################Thresholds
pes<-cbind(pm[164:167], 1, hpd[164:167,1], 2, hpd[164:167,2], 4)
rownames(pes)<-c("gamma_1", "gamma_2", "gamma_3", "gamma_4")
stargazer(pes, digits=2)

##########################################################
################Item NA analyses, Appendix K##############
##########################################################

########################################Experimental
####Function to estimate probability respondent provides experimental NA across MCMC draws with 
#different ethnic and linguistic characteristics, holding everything else constant at mean or mode
bpred<-function(cont,t1,t2,t3,beta,inc,age){
  mu<-matrix(nrow=nrow(cont),ncol=36)
  pred<-matrix(nrow=nrow(cont),ncol=36)
  for(i in 1:nrow(cont)){
    ###monolingual Russians
    mu[i,1] <-cont[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,2] <-t1[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,3] <-t2[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,4] <-t3[i,]%*%rbind(1,0,0,0,0,0)
    ###bilingual R/M Russians
    mu[i,5] <-cont[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,6] <-t1[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,7] <-t2[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,8] <-t3[i,]%*%rbind(1,0,1,0,0,0)
    ###bilingual R/U Russians
    mu[i,9] <-cont[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,10] <-t1[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,11] <-t2[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,12] <-t3[i,]%*%rbind(1,0,0,1,0,0)
    ###monolingual R Moldovans
    mu[i,13] <-cont[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,14] <-t1[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,15] <-t2[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,16] <-t3[i,]%*%rbind(1,0,0,0,1,0)
    ###bilingual R/M Moldovans
    mu[i,17] <-cont[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,18] <-t1[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,19] <-t2[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,20] <-t3[i,]%*%rbind(1,0,1,0,1,0)
    ###monolingual M Moldovans 
    mu[i,21] <-cont[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,22] <-t1[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,23] <-t2[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,24] <-t3[i,]%*%rbind(1,1,1,0,1,0)
    ###monolingual R Ukrainians
    mu[i,25] <-cont[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,26] <-t1[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,27] <-t2[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,28] <-t3[i,]%*%rbind(1,0,0,0,1,1)
    ###bilingual R/U Ukrainians
    mu[i,29] <-cont[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,30] <-t1[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,31] <-t2[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,32] <-t3[i,]%*%rbind(1,0,0,1,1,1)
    ###monolingual U Ukrainians
    mu[i,33] <-cont[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,34] <-t1[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,35] <-t2[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,36] <-t3[i,]%*%rbind(1,1,0,1,1,1)
    for(j in 1:36){
      pred[i,j]<-  pnorm(mu[i,j] + beta[i,]%*%rbind(inc[i],age))
    }
  }
  return(pred)
}


###############Coefficients for different conditions
cont<- as.matrix(cds[,seq(168,191,4)])
t1<- as.matrix(cds[,seq(169,191,4)])
t2<- as.matrix(cds[,seq(170,191,4)])
t3<- as.matrix(cds[,seq(171,191,4)])
###############Treatment-invariant coefficients
beta<- as.matrix(cds[,196:197])

############Estimate probabilities
mu<-bpred(cont,t1,t2,t3,beta,inc,age)

####Create dataset
mum<-apply(mu,2,median)
muhpd<-HPDinterval(as.mcmc(mu),prob=.9)

Condition<-rep(c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian"),9)
Ethnicity<-c(rep("Russian",12),rep("Moldovan",12),rep("Ukrainian",12))
Language<-c(rep("Russian monolingual",4),rep("Russian/Moldovan bilingual",4),
            rep("Russian/Ukrainian bilingual",4), rep("Russian monolingual",4),
            rep("Russian/Moldovan bilingual",4), rep("Moldovan monolingual",4),
            rep("Russian monolingual",4),rep("Russian/Ukrainian bilingual",4),
            rep("Ukrainian monolingual",4))


mu<-data.frame(mum,muhpd,Condition, Ethnicity,Language)
names(mu)[1]<-"value"

mu$levels<-interaction(mu$Ethnicity,mu$Language,mu$Condition)
mu$levels<-droplevels(mu$levels)
mu$levels = factor(mu$levels,sort(levels(mu$levels)))

mu$levels = factor(mu$levels,levels(mu$levels)[c(13:16,21:24,17:20,25:36,5:12,1:4)])

levels(mu$levels)<-c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian", 
                     " Control", " T_Moldovan", " T_Russian", " T_Ukrainian",
                     "  Control", "  T_Moldovan", "  T_Russian", "  T_Ukrainian",
                     "   Control", "   T_Moldovan", "   T_Russian", "   T_Ukrainian",
                     "    Control", "    T_Moldovan", "    T_Russian", "    T_Ukrainian",
                     "     Control", "     T_Moldovan", "     T_Russian", "     T_Ukrainian",
                     "      Control", "      T_Moldovan", "      T_Russian", "      T_Ukrainian",
                     "       Control", "       T_Moldovan", "       T_Russian", "       T_Ukrainian",
                     "        Control", "        T_Moldovan", "        T_Russian", "        T_Ukrainian")

mu$Ethnicity<-as.character(mu$Ethnicity)

######################
######Figure K.3######
######################

ggplot(data = mu[which(mu$Ethnicity=="Moldovan"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

ggplot(data = mu[which(mu$Ethnicity=="Ukrainian"),], 
       aes(levels,value,shape=Language,  color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="", linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

ggplot(data = mu[which(mu$Ethnicity=="Russian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="") + guides(colour = guide_legend(nrow = 2))+ 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

#####################Non-experimental NA
####Function to estimate probability respondent provides NA for given outcome across MCMC draws 
#with different ethnic and linguistic characteristics, holding everything else constant at mean
#or mode
bpreds<-function(cont,inc,age){
  mu<-matrix(nrow=nrow(cont),ncol=9)
  pred<-matrix(nrow=nrow(cont),ncol=9)
  for(i in 1:nrow(cont)){
    ###monolingual Russians
    mu[i,1] <-cont[i,]%*%rbind(1,0,0,0,0,0,inc[i],age)
    ###bilingual R/M Russians
    mu[i,2] <-cont[i,]%*%rbind(1,0,1,0,0,0,inc[i],age)
    ###bilingual R/U Russians
    mu[i,3] <-cont[i,]%*%rbind(1,0,0,1,0,0,inc[i],age)
    ###monolingual R Moldovans
    mu[i,4] <-cont[i,]%*%rbind(1,0,0,0,1,0,inc[i],age)
    ###bilingual R/M Moldovans
    mu[i,5] <-cont[i,]%*%rbind(1,0,1,0,1,0,inc[i],age)
    ###monolingual M Moldovans 
    mu[i,6] <-cont[i,]%*%rbind(1,1,1,0,1,0,inc[i],age)
    ###monolingual R Ukrainians
    mu[i,7] <-cont[i,]%*%rbind(1,0,0,0,1,1,inc[i],age)
    ###bilingual R/U Ukrainians
    mu[i,8] <-cont[i,]%*%rbind(1,0,0,1,1,1,inc[i],age)
    ###monolingual U Ukrainians
    mu[i,9] <-cont[i,]%*%rbind(1,1,0,1,1,1,inc[i],age)
    for(j in 1:9){
      pred[i,j]<-  pnorm(mu[i,j])
    }
  }
  return(pred)
}


########################Different outcome coefficient estimates
####################Integration w Moldova 
contMol<- as.matrix(cds[,c(seq(208,267,10),seq(278,297,10))])
####################Autonomy in Moldova 
contAut<- as.matrix(cds[,c(seq(209,267,10),seq(279,297,10))])
####################Confederation w Moldova 
contCon<- as.matrix(cds[,c(seq(210,267,10),seq(280,297,10))])
##################Independence 
contInd<- as.matrix(cds[,c(seq(211,267,10),seq(281,297,10))])
#################Integration w Russia 
contRus<- as.matrix(cds[,c(seq(212,267,10),seq(282,297,10))])
######################Autonomy in Russia 
contRA<- as.matrix(cds[,c(seq(213,267,10),seq(283,297,10))])
#################Confederation w Russia 
contRC<- as.matrix(cds[,c(seq(214,267,10),seq(284,297,10))])
#################Integration w Ukraine
contUkr<- as.matrix(cds[,c(seq(215,267,10),seq(285,297,10))])
######################Autonomy in Ukraine 
contUA<- as.matrix(cds[,c(seq(216,267,10),seq(286,297,10))])
#################Confederation w Ukraine
contUC<- as.matrix(cds[,c(seq(217,267,10),seq(287,297,10))])

#################Figure K.1
############To create graphic for relevant outcome, insert outcome matrix here
mu<-bpreds(contInd,inc,age)
mum<-apply(mu,2,median)
muhpd<-HPDinterval(as.mcmc(mu),prob=.9)


Language<-as.factor(c("Russian monolingual","Russian/Moldovan bilingual","Russian/Ukrainian bilingual",
                      "Russian monolingual", "Russian/Moldovan bilingual", "Moldovan monolingual",
                      "Russian monolingual", "Russian/Ukrainian bilingual", "Ukrainian monolingual"))

Ethnicity <- as.factor(c(rep("Russian", 3),rep("Moldovan", 3),rep("Ukrainian", 3)))

mu<-data.frame(mum,muhpd,Language,Ethnicity)

names(mu)[1]<-"value"

mu$levels<-interaction(mu$Ethnicity,mu$Language)
mu$levels<-droplevels(mu$levels)
mu$levels = factor(mu$levels,levels(mu$levels)[c(3,7,6,4,8,9,2,5,1)])
levels(mu$levels)<-c("  Russian monolingual", "  Russian/Ukrainian bilingual", 
                     "  Russian/Moldovan bilingual", " Russian monolingual",        
                     " Russian/Ukrainian bilingual", " Ukrainian monolingual",      
                     "Russian monolingual", "Russian/Moldovan bilingual",  
                     "Moldovan monolingual")

mu$Ethnicity<-as.character(mu$Ethnicity)


ggplot(data = mu, aes(levels,value,color=Ethnicity, shape=Ethnicity))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(3.5,6.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="") + 
  theme(text = element_text(size=15), legend.position=c(.75,.25))  + 
  coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

##################################################
######################Language salience statistics
#############Table B.1############################
##################################################

##############PMR government support for language
###Russian
mean(pmrdata$pgs_r>3, na.rm=T)
###Gagauz
mean(pmrdata$pgs_u>3, na.rm=T)
###Moldovan
mean(pmrdata$pgs_m>3, na.rm=T)

##############Moldovan government support for language
###Russian
mean(pmrdata$mgs_r>3, na.rm=T)
###Ukrainian
mean(pmrdata$mgs_u>3, na.rm=T)
###Moldovan
mean(pmrdata$mgs_m>3, na.rm=T)

####################Gov comparison in txt
mean(pmrdata$pgs_r>pmrdata$mgs_r,na.rm=T)
mean(pmrdata$pgs_u>pmrdata$mgs_u,na.rm=T)
mean(pmrdata$pgs_m>pmrdata$mgs_m,na.rm=T)
mean(pmrdata$pgs_m<pmrdata$mgs_m,na.rm=T)

##############Frequency of language use in PMR
###Russian
mean(pmrdata$plf_r>3, na.rm=T)
###Ukrainian
mean(pmrdata$plf_u>3, na.rm=T)
###Moldovan
mean(pmrdata$plf_m>3, na.rm=T)

##############Frequency of language use in Moldova
###Russian
mean(pmrdata$mlf_r>3, na.rm=T)
###Ukrainian
mean(pmrdata$mlf_u>3, na.rm=T)
###Moldovan
mean(pmrdata$mlf_m>3, na.rm=T)

####################Frequency of use comparison in txt
mean(pmrdata$plf_r>pmrdata$mlf_r,na.rm=T)
mean(pmrdata$plf_u>pmrdata$mlf_u,na.rm=T)
mean(pmrdata$plf_m<pmrdata$mlf_m,na.rm=T)



#################################################################
##############################Components of ethnic identity stats
###########################Table C1##############################

###########Create data frame w components by ethnic group
###Note-NAs/DKs coded as 0
eth_id<-data.frame(Component=c("Mother's nationality", "Father's nationality", 
                               "Following cultural traditions", "Knowledge of language", 
                               "Religious belief", "Physical appearance", "Self-identification", 
                               "Location of residence", "Citizenship", "Common history",
                               "Professional traditions", "Marrying co-ethnic", 
                               "Educating children in language"),
                   Russian=c(mean(pmrdata$N3a[pmrdata$Russian==1]), 
                             mean(pmrdata$N3b[pmrdata$Russian==1]),
                             mean(pmrdata$N3c[pmrdata$Russian==1]),
                             mean(pmrdata$N3d[pmrdata$Russian==1]),
                             mean(pmrdata$N3e[pmrdata$Russian==1]),
                             mean(pmrdata$N3f[pmrdata$Russian==1]),
                             mean(pmrdata$N3g[pmrdata$Russian==1]),
                             mean(pmrdata$N3h[pmrdata$Russian==1]),
                             mean(pmrdata$N3i[pmrdata$Russian==1]),
                             mean(pmrdata$N3j[pmrdata$Russian==1]),
                             mean(pmrdata$N3k[pmrdata$Russian==1]),
                             mean(pmrdata$N3l[pmrdata$Russian==1]),
                             mean(pmrdata$N3m[pmrdata$Russian==1])),
                   Moldovan=c(mean(pmrdata$N3a[pmrdata$Moldovan==1]), 
                              mean(pmrdata$N3b[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3c[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3d[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3e[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3f[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3g[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3h[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3i[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3j[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3k[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3l[pmrdata$Moldovan==1]),
                              mean(pmrdata$N3m[pmrdata$Moldovan==1])),
                   Ukrainian=c(mean(pmrdata$N3a[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3b[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3c[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3d[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3e[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3f[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3g[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3h[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3i[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3j[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3k[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3l[pmrdata$Ukrainian==1]),
                               mean(pmrdata$N3m[pmrdata$Ukrainian==1])),
                   Overall=c(mean(pmrdata$N3a), 
                             mean(pmrdata$N3b),
                             mean(pmrdata$N3c),
                             mean(pmrdata$N3d),
                             mean(pmrdata$N3e),
                             mean(pmrdata$N3f),
                             mean(pmrdata$N3g),
                             mean(pmrdata$N3h),
                             mean(pmrdata$N3i),
                             mean(pmrdata$N3j),
                             mean(pmrdata$N3k),
                             mean(pmrdata$N3l),
                             mean(pmrdata$N3m)))

####Table, ordered by overall order
eth_id[order(eth_id$Overall, decreasing=T),]


#######################################
#########Age of acquisition, Appendix D
#######################################

############Figure D.1

##########Russian
qplot(aa_r,1-rf,data=pmrdata, xlab="Age of acquisition", ylab="Fluent in Russian") + 
  geom_jitter(height=.01) + stat_smooth(method="lm") + theme_bw() + 
  coord_cartesian(ylim=c(0,1),xlim=c(0,50))

###########Moldovan
qplot(aa_m,mf,data=pmrdata, xlab="Age of acquisition", ylab="Fluent in Moldovan") + 
  geom_jitter(height=.01) + stat_smooth(method="lm") + theme_bw() + 
  coord_cartesian(ylim=c(0,1),xlim=c(0,50))

##########Ukrainian
qplot(aa_u,uf,data=pmrdata, xlab="Age of acquisition", ylab="Fluent in Ukrainian") + 
  geom_jitter(height=.01) + stat_smooth(method="lm") + theme_bw() + 
  coord_cartesian(ylim=c(0,1),xlim=c(0,50))

################################################
#######Reasons for learning language, Appendix E
################################################

###############Russian
#####Reason for learning language x fluent
raa<-table(pmrdata$reason_r,1-pmrdata$rf)
raa[,2]/apply(raa,1,sum) ###Proportion fluent

####Family Reason for learning language x fluent
raa<-table(pmrdata$freason_r,1-pmrdata$rf)
raa[,2]/apply(raa,1,sum) ###Proportion fluent

###############Moldovan
#####Reason for learning language x fluent
maa<-table(pmrdata$reason_m,pmrdata$mf)
maa[,2]/apply(maa,1,sum) ###Proportion fluent

####Family Reason for learning language x fluent
maa<-table(pmrdata$freason_m,pmrdata$mf)
maa[,2]/apply(maa,1,sum) ###Proportion fluent

###############Ukrainian
#####Reason for learning language x fluent
uaa<-table(pmrdata$reason_u,pmrdata$uf)
uaa[,2]/apply(uaa,1,sum) ###Proportion fluent

####Family Reason for learning language x fluent
uaa<-table(pmrdata$freason_u,pmrdata$uf)
uaa[,2]/apply(uaa,1,sum) ###Proportion fluent

############################################
##########Descriptive statistics, Appendix G
############################################

###Table G.1
#####Language
mean(pmrdata$rf,na.rm=T)
sum(!is.na(pmrdata$rf))

mean(pmrdata$mf,na.rm=T)
sum(!is.na(pmrdata$mf))

mean(pmrdata$uf,na.rm=T)
sum(!is.na(pmrdata$uf))

####Ethnicity (no NA)
mean(pmrdata$Russian)
mean(pmrdata$Moldovan)
mean(pmrdata$Ukrainian)
mean(pmrdata$Other)

####Income
income<- apply(as.matrix(cds[,396:972]),2,mean) ####Estimate average income across draws for respondents
mean(income)
sd(income)

###Age
mean(pmrdata$age)
sd(pmrdata$age)

###Higher education
mean(pmrdata$highered,na.rm=T)
sum(!is.na(pmrdata$highered))

##Male 
mean(pmrdata$male)

##Urban
mean(pmrdata$urban)

#########Table G.2

###Independence
median(pmrdata$S2d,na.rm=T)
quantile(pmrdata$S2d,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$S2d))
sum(pmrdata$S2d_dk)/sum(is.na(pmrdata$S2d))

###Integration with Moldova
median(pmrdata$S2a,na.rm=T)
quantile(pmrdata$S2a,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$S2a))
sum(pmrdata$S2a_dk)/sum(is.na(pmrdata$S2a))

median(pmrdata$S2b,na.rm=T)
quantile(pmrdata$S2b,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$S2b))
sum(pmrdata$S2b_dk)/sum(is.na(pmrdata$S2b))

median(pmrdata$S2c,na.rm=T)
quantile(pmrdata$S2c,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$S2c))
sum(pmrdata$S2c_dk)/sum(is.na(pmrdata$S2c))

####One Moldova response
577- sum(is.na(pmrdata$S2a) & is.na(pmrdata$S2b) & is.na(pmrdata$S2c))

###Integration with Russia
median(pmrdata$R1a,na.rm=T)
quantile(pmrdata$R1a,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$R1a))
sum(pmrdata$R1a_dk)/sum(is.na(pmrdata$R1a))

################Autonomy w/in Russia
median(pmrdata$R1b,na.rm=T)
quantile(pmrdata$R1b,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$R1b))
sum(pmrdata$R1b_dk)/sum(is.na(pmrdata$R1b))

############Confederation w Russia
median(pmrdata$R1c,na.rm=T)
quantile(pmrdata$R1c,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$R1c))
sum(pmrdata$R1c_dk)/sum(is.na(pmrdata$R1c))

###One Russia response
577-sum(is.na(pmrdata$R1a) & is.na(pmrdata$R1b) & is.na(pmrdata$R1c))

###Integration with Ukraine
median(pmrdata$U1a,na.rm=T)
quantile(pmrdata$U1a,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$U1a))
sum(pmrdata$U1a_dk)/sum(is.na(pmrdata$U1a))

################Autonomy w/in Ukraine
median(pmrdata$U1b,na.rm=T)
quantile(pmrdata$U1b,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$U1b))
sum(pmrdata$U1b_dk)/sum(is.na(pmrdata$U1b))

############Confederation w Ukraine
median(pmrdata$U1c,na.rm=T)
quantile(pmrdata$U1c,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$U1c))
sum(pmrdata$U1c_dk)/sum(is.na(pmrdata$U1c))

###One Russia response
577-sum(is.na(pmrdata$U1a) & is.na(pmrdata$U1b) & is.na(pmrdata$U1c))

#####1 Likert scale complete
577-sum(is.na(pmrdata$S2d) & is.na(pmrdata$S2a) & is.na(pmrdata$S2b) & is.na(pmrdata$S2c) & is.na(pmrdata$R1a) & is.na(pmrdata$R1b) & is.na(pmrdata$R1c) & is.na(pmrdata$U1a) & is.na(pmrdata$U1b) & is.na(pmrdata$U1c))
#####1 outcome
577-sum(is.na(pmrdata$DV1) & is.na(pmrdata$S2d) & is.na(pmrdata$S2a) & is.na(pmrdata$S2b) & is.na(pmrdata$S2c) & is.na(pmrdata$R1a) & is.na(pmrdata$R1b) & is.na(pmrdata$R1c) & is.na(pmrdata$U1a) & is.na(pmrdata$U1b) & is.na(pmrdata$U1c))


###Table G.3
###########Outcome
median(pmrdata$DV1,na.rm=T)
quantile(pmrdata$DV1,na.rm=T)[c(2,4)]
sum(!is.na(pmrdata$DV1))
sum(pmrdata$DV1_dk)/sum(is.na(pmrdata$DV1))

###Treatments
table(pmrdata$tg)/577


############################################
#########################T-tests, Appendix J
############################################

###########Table J.1
##########Overall
t.test(pmrdata$DV1[which(pmrdata$tg=="Control")] > 3, pmrdata$DV1[which(pmrdata$tg=="T1")] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control")] > 3, pmrdata$DV1[which(pmrdata$tg=="T2")] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control")] > 3, pmrdata$DV1[which(pmrdata$tg=="T3")] > 3)

####Fluent in Russian
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$rf==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$rf==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$rf==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$rf==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$rf==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$rf==0)] > 3)

####Not fluent in Russian
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$rf==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$rf==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$rf==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$rf==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$rf==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$rf==1)] > 3)

###Fluent in Moldovan
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$mf==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$mf==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$mf==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$mf==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$mf==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$mf==1)] > 3)

###Not fluent in Moldovan
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$mf==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$mf==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$mf==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$mf==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$mf==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$mf==0)] > 3)

###Fluent in Ukrainian
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$uf==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$uf==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$uf==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$uf==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$uf==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$uf==1)] > 3)

###Not fluent in Ukrainian
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$uf==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$uf==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$uf==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$uf==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$uf==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$uf==0)] > 3)

####Russian
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Russian==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$Russian==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Russian==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$Russian==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Russian==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$Russian==1)] > 3)

###Not Russian
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Russian==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$Russian==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Russian==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$Russian==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Russian==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$Russian==0)] > 3)

###Ethnic Moldovan
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Moldovan==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$Moldovan==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Moldovan==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$Moldovan==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Moldovan==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$Moldovan==1)] > 3)

#####Not Moldovan
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Moldovan==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$Moldovan==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Moldovan==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$Moldovan==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Moldovan==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$Moldovan==0)] > 3)

####Ethnic Ukrainian
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Ukrainian==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$Ukrainian==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Ukrainian==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$Ukrainian==1)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Ukrainian==1)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$Ukrainian==1)] > 3)

###Not Ukrainian
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Ukrainian==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T1" & pmrdata$Ukrainian==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Ukrainian==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T2" & pmrdata$Ukrainian==0)] > 3)
t.test(pmrdata$DV1[which(pmrdata$tg=="Control" & pmrdata$Ukrainian==0)] > 3, pmrdata$DV1[which(pmrdata$tg=="T3" & pmrdata$Ukrainian==0)] > 3)


###########################################################
#####Analysis with random experimental clustering variation
###############################################Appendix J.6
###########################################################

######################Experimental threshold initial values
taustar<-polr(method="probit", as.factor(DV1) ~ rf + mf + uf + NRussian + Ukrainian + Other + income + 
                age + male + highered + urban, data=pmrdata)$zeta

###########Data file
hmod<-list(N=nrow(pmrdata), tg=pmrdata$tg, 
           y1=pmrdata$DV1, y2=as.matrix(pmrdata[,24:33]), 
           y1NA=pmrdata$DV1na, y2NA=as.matrix(pmrdata[,34:43]), 
           g0=rep(0,4), G0=diag(1,4), 
           t0=rep(0,3), T0=diag(.01,3), 
           rf=pmrdata$rf, mf=pmrdata$mf, uf=pmrdata$uf, NRussian=pmrdata$NRussian, 
           Other=pmrdata$Other, Ukrainian=pmrdata$Ukrainian, Male=pmrdata$male, 
           Age=pmrdata$age, Urban=pmrdata$urban,
           NIOTA=12, I0=diag(1,12), i0=rep(0,12), 
           NALPHA=7, A0=diag(1,7), a0=rep(0,7), 
           B0=diag(1,5), b0=rep(0,5), 
           HigherEd=pmrdata$highered, Married=pmrdata$married, Locality=pmrdata$locality, 
           Prof=pmrdata$prof, Resp=pmrdata$resp, House=pmrdata$house, ymin=pmrdata$ymin, 
           ymax=pmrdata$ymax, 
           bI0=rep(0,5), BI0=diag(1,5), ae1=rep(1,14))

#hmodel <- run.jags(method="parallel",model="pF2_exp.txt", inits=list(taustar=taustar), 
#                   monitor=c("alpha", "beta", "tau", "mu.a",
#                             "tauExp", 
#                             "pr", "pm",  "pu",
#                             "bI", "betaI","gammaI","alphaI",
#                             "pIe", "pIm","e1", "rIl","rIh",
#                             "Income"), data=hmod, n.chains=4, adapt=5000, burnin=5000000, 
#                   sample=1000, thin=50000, summarise=F,plots=F, modules=c("glm", "lecuyer"))

#save(hmodel, file = "pF1_exp.RData")
load("pF1_exp.RData")

###Check convergence,  "Other" ethnic identity due to small sample size across experimental conditions
h1<-as.mcmc.list(hmodel)
gd<-gelman.diag(h1[,1:118],multivariate = F)
sum(gd[[1]][,1]>1.1)

####Convert to single chain
cds<-as.mcmc(hmodel)

###Mean income and age
inc<- apply(as.matrix(cds[,119:695]),1,mean)
age<-mean(pmrdata$age)


############################experimental analysis
####Function to estimate probability respondent supports PMR independence across MCMC draws 
#with different ethnic and linguistic characteristics, holding everything else constant at mean
#or mode
bpred<-function(cont,t1,t2,t3,eps,beta,inc,age){
  mu<-matrix(nrow=nrow(cont),ncol=36)
  pred<-matrix(nrow=nrow(cont),ncol=36)
  for(i in 1:nrow(cont)){
    ###monolingual Russians
    mu[i,1] <-cont[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,2] <-t1[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,3] <-t2[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,4] <-t3[i,]%*%rbind(1,0,0,0,0,0)
    ###bilingual R/M Russians
    mu[i,5] <-cont[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,6] <-t1[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,7] <-t2[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,8] <-t3[i,]%*%rbind(1,0,1,0,0,0)
    ###bilingual R/U Russians
    mu[i,9] <-cont[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,10] <-t1[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,11] <-t2[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,12] <-t3[i,]%*%rbind(1,0,0,1,0,0)
    ###monolingual R Moldovans
    mu[i,13] <-cont[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,14] <-t1[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,15] <-t2[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,16] <-t3[i,]%*%rbind(1,0,0,0,1,0)
    ###bilingual R/M Moldovans
    mu[i,17] <-cont[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,18] <-t1[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,19] <-t2[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,20] <-t3[i,]%*%rbind(1,0,1,0,1,0)
    ###monolingual M Moldovans 
    mu[i,21] <-cont[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,22] <-t1[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,23] <-t2[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,24] <-t3[i,]%*%rbind(1,1,1,0,1,0)
    ###monolingual R Ukrainians
    mu[i,25] <-cont[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,26] <-t1[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,27] <-t2[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,28] <-t3[i,]%*%rbind(1,0,0,0,1,1)
    ###bilingual R/U Ukrainians
    mu[i,29] <-cont[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,30] <-t1[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,31] <-t2[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,32] <-t3[i,]%*%rbind(1,0,0,1,1,1)
    ###monolingual U Ukrainians
    mu[i,33] <-cont[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,34] <-t1[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,35] <-t2[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,36] <-t3[i,]%*%rbind(1,1,0,1,1,1)
    for(j in 1:36){
      pred[i,j]<-  1- pnorm(eps[i]  - (mu[i,j] + beta[i,]%*%rbind(inc[i],age)))
    }
  }
  return(pred)
}

####Third threshold
eps<-as.matrix(cds[,36])

###Experimental coefficients
cont<- as.matrix(cds[,seq(1,24,4)])
t1<- as.matrix(cds[,seq(2,24,4)])
t2<- as.matrix(cds[,seq(3,24,4)])
t3<- as.matrix(cds[,seq(4,24,4)])
beta<- as.matrix(cds[,29:30])

###Treatment-invariant coefficients
mu<-bpred(cont,t1,t2,t3,eps,beta,inc,age)

###Median and HPD interval
mum<-apply(mu,2,median)
muhpd<-HPDinterval(as.mcmc(mu),prob=.95)

########Create dataset
Condition<-rep(c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian"),9)
Ethnicity<-c(rep("Russian",12),rep("Moldovan",12),rep("Ukrainian",12))
Language<-factor(c(rep("Russian monolingual",4),rep("Russian/Moldovan bilingual",4),
            rep("Russian/Ukrainian bilingual",4), rep("Russian monolingual",4),
            rep("Russian/Moldovan bilingual",4), rep("Moldovan monolingual",4),
            rep("Russian monolingual",4),rep("Russian/Ukrainian bilingual",4),
            rep("Ukrainian monolingual",4)))

mu<-data.frame(mum,muhpd,Condition, Ethnicity,Language)
names(mu)[1]<-"value"

###Subgroup
mu$levels<-interaction(mu$Ethnicity,mu$Language,mu$Condition)
mu$levels<-droplevels(mu$levels)
mu$levels = factor(mu$levels,sort(levels(mu$levels)))

mu$levels = factor(mu$levels,levels(mu$levels)[c(13:16,21:24,17:20,25:36,5:12,1:4)])

levels(mu$levels)<-c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian", 
                     " Control", " T_Moldovan", " T_Russian", " T_Ukrainian",
                     "  Control", "  T_Moldovan", "  T_Russian", "  T_Ukrainian",
                     "   Control", "   T_Moldovan", "   T_Russian", "   T_Ukrainian",
                     "    Control", "    T_Moldovan", "    T_Russian", "    T_Ukrainian",
                     "     Control", "     T_Moldovan", "     T_Russian", "     T_Ukrainian",
                     "      Control", "      T_Moldovan", "      T_Russian", "      T_Ukrainian",
                     "       Control", "       T_Moldovan", "       T_Russian", "       T_Ukrainian",
                     "        Control", "        T_Moldovan", "        T_Russian", "        T_Ukrainian")

mu$Ethnicity<-as.character(mu$Ethnicity)

####Figure J.7
ggplot(data = mu[which(mu$Ethnicity=="Russian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

###Figure J.8
ggplot(data = mu[which(mu$Ethnicity=="Moldovan"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

##Figure J.9
ggplot(data = mu[which(mu$Ethnicity=="Ukrainian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))


############################################
###############Russian influence, Appendix L
############################################

###Proportion supporting integration with Russia w/o autonomy
mean(pmrdata$R1a>3,na.rm=T)

###Proportion supporting independence
mean(pmrdata$S2d>3,na.rm=T)

####Figure L.1
qplot(log(Russia+1),1-rf, xlab="Time in Russia", ylab="Fluent in Russian", data=pmrdata) + 
  stat_smooth(method="lm") + theme_bw() +  coord_cartesian(ylim=c(-.1,1.1)) + scale_color_grey() + 
  theme(text = element_text(size=15), legend.position  = "top") + geom_jitter(height=.05,width=.05)

####Figure L.2
qplot(log(Russia+1),S2d, xlab="Time in Russia", ylab="Support for independence", data=pmrdata) + 
  stat_smooth(method="lm") + scale_color_grey() + theme_bw() + 
  theme(text = element_text(size=15), legend.position = "top") +  geom_jitter(height=.1,width=.05) + 
  coord_cartesian(ylim=c(1,5)) 

qplot(log(Russia+1),R1a, xlab="Time in Russia", ylab="Support for integration with Russia", data=pmrdata) + 
  stat_smooth(method="lm") + scale_color_grey() + theme_bw() + 
  theme(text = element_text(size=15), legend.position = "top") +  geom_jitter(height=.1,width=.05) + 
  coord_cartesian(ylim=c(1,5)) 

##########################################
##Media, Appendix M#######################
##########################################

##############Table M.1
table(1-pmrdata$rf, 1-pmrdata$rf_m)/sum(!is.na(pmrdata$rf) & !is.na(pmrdata$rf_m) )
table(pmrdata$mf, pmrdata$mf_m)/sum(!is.na(pmrdata$mf) & !is.na(pmrdata$mf_m) )
table(pmrdata$uf, pmrdata$uf_m)/sum(!is.na(pmrdata$uf) & !is.na(pmrdata$uf_m) )

#####Figure M.2
qplot(as.factor(1-rf), S2a, color=as.factor(mf_m), fill=as.factor(mf_m), data=pmrdata[which(!is.na(pmrdata$rf) & !is.na(pmrdata$mf_m)),],
      xlab="Fluent in Russian", ylab="Supports integration with Moldova", geom="boxplot") + labs(color="Watches TV in Moldovan", fill="Watches TV in Moldovan") + 
  theme_bw() + scale_color_grey() + scale_fill_grey() + theme(text = element_text(size=15), legend.position = "top")

qplot(as.factor(1-rf), S2d, color=as.factor(mf_m), fill=as.factor(mf_m), data=pmrdata[which(!is.na(pmrdata$rf) & !is.na(pmrdata$mf_m)),],
      xlab="Fluent in Russian", ylab="Supports independence", geom="boxplot") + labs(color="Watches TV in Moldovan", fill="Watches TV in Moldovan") + 
  theme_bw() + scale_color_grey() + scale_fill_grey() + theme(text = element_text(size=15), legend.position = "top")

qplot(as.factor(1-rf), S2d, color=as.factor(uf_m), fill=as.factor(uf_m), data=pmrdata[which(!is.na(pmrdata$rf) & !is.na(pmrdata$uf_m)),],
      xlab="Fluent in Russian", ylab="Supports independence", geom="boxplot") + labs(color="Watches TV in Ukrainian", fill="Watches TV in Ukrainian") + 
  theme_bw() + scale_color_grey() + scale_fill_grey() + theme(text = element_text(size=15), legend.position = "top")

qplot(as.factor(1-rf), R1a, color=as.factor(mf_m), fill=as.factor(mf_m), data=pmrdata[which(!is.na(pmrdata$rf) & !is.na(pmrdata$mf_m)),],
      xlab="Fluent in Russian", ylab="Supports integration with Russia", geom="boxplot") + labs(color="Watches TV in Moldovan", fill="Watches TV in Moldovan") + 
  theme_bw() + scale_color_grey() + scale_fill_grey() + theme(text = element_text(size=15), legend.position = "top")

qplot(as.factor(1-rf), R1a, color=as.factor(uf_m), fill=as.factor(uf_m), data=pmrdata[which(!is.na(pmrdata$rf) & !is.na(pmrdata$uf_m)),],
      xlab="Fluent in Russian", ylab="Supports integration with Russia", geom="boxplot") + labs(color="Watches TV in Ukrainian", fill="Watches TV in Ukrainian") + 
  theme_bw() + scale_color_grey() + scale_fill_grey() + theme(text = element_text(size=15), legend.position = "top")

################################################################
##########Regression analyses of media on support for separatism
################################################################

##############Experimental threshold initial values
taustar<-polr(method="probit", as.factor(DV1) ~ rf_m + mf_m + uf_m + NRussian + Ukrainian + Other + income + 
                age + male + highered + urban, data=pmrdata)$zeta


##############Observational threshold initial values
gammastar<-matrix(nrow=10,ncol=4)
for(i in 1:10){
  gammastar[i,]<-polr(method="probit", as.factor(pmrdata[,i+23]) ~ pmrdata$rf_m + pmrdata$mf_m + 
                        pmrdata$uf_m + pmrdata$NRussian + pmrdata$Ukrainian + pmrdata$Other + 
                        pmrdata$income + pmrdata$age + pmrdata$male + pmrdata$highered + 
                        pmrdata$urban)$zeta
}
gammastar<-apply(gammastar,2,mean)

###########Data file
hmod<-list(N=nrow(pmrdata), tg=pmrdata$tg, 
           y1=pmrdata$DV1, y2=as.matrix(pmrdata[,24:33]), 
           g0=rep(0,4), G0=diag(.01,4), 
           t0=rep(0,3), T0=diag(.01,3), 
           rf=pmrdata$rf_m, mf=pmrdata$mf_m, uf=pmrdata$uf_m, NRussian=pmrdata$NRussian, 
           Other=pmrdata$Other, Ukrainian=pmrdata$Ukrainian, Male=pmrdata$male, 
           Age=pmrdata$age, Urban=pmrdata$urban,
           NIOTA=12, I0=diag(1,12), i0=rep(0,12), 
           NALPHA=7, A0=diag(1,7), a0=rep(0,7), 
           B0=diag(1,5), b0=rep(0,5), 
           HigherEd=pmrdata$highered, Married=pmrdata$married, Locality=pmrdata$locality, 
           Prof=pmrdata$prof, Resp=pmrdata$resp, House=pmrdata$house, ymin=pmrdata$ymin, 
           ymax=pmrdata$ymax, 
           bI0=rep(0,5), BI0=diag(1,5), ae1=rep(1,14))

#hmodel <- run.jags(method="parallel",model="pFr.txt", inits=list(gammastar=gammastar,
#                                                                 taustar=taustar), 
#                   monitor=c("alpha", "beta", "tau", "mu.a",
#                             "iota", "gamma",
#                             "pr", "pm",  "pu",
#                             "bI", "betaI","gammaI","alphaI",
#                             "pIe", "pIm","e1", "rIl","rIh",
#                             "Income"), data=hmod, n.chains=4, adapt=5000, burnin=50000, 
#                   sample=1000, thin=500, summarise=F,plots=F, modules=c("glm", "lecuyer"))
#save(hmodel, file = "pF1r.RData")
load("pF1r.RData")

#####Check convergence
h1<-as.mcmc.list(hmodel)
gd<-gelman.diag(h1[,1:235],multivariate=F)
sum(gd[[1]][,1]>1.1)

#####Convert to single chain
cds<-as.mcmc(hmodel)

####Mean income and age across posterior draws
inc<- apply(as.matrix(cds[,236:812]),1,mean)
age<-mean(pmrdata$age)

############################experimental analysis
##############Not reported in text
####Function to estimate probability respondent supports PMR independence across MCMC draws 
#with different ethnic and linguistic characteristics, holding everything else constant at mean
#or mode
bpred<-function(cont,t1,t2,t3,eps,beta,inc,age){
  mu<-matrix(nrow=nrow(cont),ncol=36)
  pred<-matrix(nrow=nrow(cont),ncol=36)
  for(i in 1:nrow(cont)){
    ###monolingual Russians
    mu[i,1] <-cont[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,2] <-t1[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,3] <-t2[i,]%*%rbind(1,0,0,0,0,0)
    mu[i,4] <-t3[i,]%*%rbind(1,0,0,0,0,0)
    ###bilingual R/M Russians
    mu[i,5] <-cont[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,6] <-t1[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,7] <-t2[i,]%*%rbind(1,0,1,0,0,0)
    mu[i,8] <-t3[i,]%*%rbind(1,0,1,0,0,0)
    ###bilingual R/U Russians
    mu[i,9] <-cont[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,10] <-t1[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,11] <-t2[i,]%*%rbind(1,0,0,1,0,0)
    mu[i,12] <-t3[i,]%*%rbind(1,0,0,1,0,0)
    ###monolingual R Moldovans
    mu[i,13] <-cont[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,14] <-t1[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,15] <-t2[i,]%*%rbind(1,0,0,0,1,0)
    mu[i,16] <-t3[i,]%*%rbind(1,0,0,0,1,0)
    ###bilingual R/M Moldovans
    mu[i,17] <-cont[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,18] <-t1[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,19] <-t2[i,]%*%rbind(1,0,1,0,1,0)
    mu[i,20] <-t3[i,]%*%rbind(1,0,1,0,1,0)
    ###monolingual M Moldovans 
    mu[i,21] <-cont[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,22] <-t1[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,23] <-t2[i,]%*%rbind(1,1,1,0,1,0)
    mu[i,24] <-t3[i,]%*%rbind(1,1,1,0,1,0)
    ###monolingual R Ukrainians
    mu[i,25] <-cont[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,26] <-t1[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,27] <-t2[i,]%*%rbind(1,0,0,0,1,1)
    mu[i,28] <-t3[i,]%*%rbind(1,0,0,0,1,1)
    ###bilingual R/U Ukrainians
    mu[i,29] <-cont[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,30] <-t1[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,31] <-t2[i,]%*%rbind(1,0,0,1,1,1)
    mu[i,32] <-t3[i,]%*%rbind(1,0,0,1,1,1)
    ###monolingual U Ukrainians
    mu[i,33] <-cont[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,34] <-t1[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,35] <-t2[i,]%*%rbind(1,1,0,1,1,1)
    mu[i,36] <-t3[i,]%*%rbind(1,1,0,1,1,1)
    for(j in 1:36){
      pred[i,j]<-  1- pnorm(eps[i]  - (mu[i,j] + beta[i,]%*%rbind(inc[i],age)))
    }
  }
  return(pred)
}


######Third threshold
eps<-as.matrix(cds[,36])

##############Experimental coefficients
cont<- as.matrix(cds[,seq(1,24,4)])
t1<- as.matrix(cds[,seq(2,24,4)])
t2<- as.matrix(cds[,seq(3,24,4)])
t3<- as.matrix(cds[,seq(4,24,4)])
###Treatment invariant coefficients
beta<- as.matrix(cds[,29:30])

#####Estimate posterior predictions
mu<-bpred(cont,t1,t2,t3,eps,beta,inc,age)

###Create dataset
mum<-apply(mu,2,median)
muhpd<-HPDinterval(as.mcmc(mu),prob=.9)

Condition<-rep(c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian"),9)
Ethnicity<-c(rep("Russian",12),rep("Moldovan",12),rep("Ukrainian",12))
Language<-factor(c(rep("Russian monolingual",4),rep("Russian/Moldovan bilingual",4),
            rep("Russian/Ukrainian bilingual",4), rep("Russian monolingual",4),
            rep("Russian/Moldovan bilingual",4), rep("Moldovan monolingual",4),
            rep("Russian monolingual",4),rep("Russian/Ukrainian bilingual",4),
            rep("Ukrainian monolingual",4)))

mu<-data.frame(mum,muhpd,Condition, Ethnicity,Language)
names(mu)[1]<-"value"

####Control
muc<-mu[mu$Condition=="Control",]
muc$levels<-interaction(muc$Ethnicity,muc$Language)
muc$levels<-droplevels(muc$levels)
muc$levels = factor(muc$levels,levels(muc$levels)[c(3,7,6,4,8,9,2,5,1)])
levels(muc$levels)<-c("  Russian monolingual", "  Russian/Ukrainian bilingual", 
                      "  Russian/Moldovan bilingual", " Russian monolingual",        
                      " Russian/Ukrainian bilingual", " Ukrainian monolingual",      
                      "Russian monolingual", "Russian/Moldovan bilingual",  
                      "Moldovan monolingual")

muc$Ethnicity<-as.character(muc$Ethnicity)

ggplot(data = muc, aes(levels,value,color=Ethnicity, shape=Ethnicity))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(3.5,6.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="") + 
  theme(text = element_text(size=15), legend.position=c(.25,.25))  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

###Subgroup
mu$levels<-interaction(mu$Ethnicity,mu$Language,mu$Condition)
mu$levels<-droplevels(mu$levels)
mu$levels = factor(mu$levels,sort(levels(mu$levels)))

mu$levels = factor(mu$levels,levels(mu$levels)[c(13:16,21:24,17:20,25:36,5:12,1:4)])

levels(mu$levels)<-c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian", 
                     " Control", " T_Moldovan", " T_Russian", " T_Ukrainian",
                     "  Control", "  T_Moldovan", "  T_Russian", "  T_Ukrainian",
                     "   Control", "   T_Moldovan", "   T_Russian", "   T_Ukrainian",
                     "    Control", "    T_Moldovan", "    T_Russian", "    T_Ukrainian",
                     "     Control", "     T_Moldovan", "     T_Russian", "     T_Ukrainian",
                     "      Control", "      T_Moldovan", "      T_Russian", "      T_Ukrainian",
                     "       Control", "       T_Moldovan", "       T_Russian", "       T_Ukrainian",
                     "        Control", "        T_Moldovan", "        T_Russian", "        T_Ukrainian")

mu$Ethnicity<-as.character(mu$Ethnicity)

ggplot(data = mu[which(mu$Ethnicity=="Russian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

ggplot(data = mu[which(mu$Ethnicity=="Ukrainian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))


ggplot(data = mu[which(mu$Ethnicity=="Moldovan"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

#######################Non-experimental analysis
#Function to estimate posterior probability respondent w different ethnic and linguistic
#characteristics supports a given separatist outcome
bpreds<-function(cont,eps,inc,age){
  mu<-matrix(nrow=nrow(cont),ncol=9)
  pred<-matrix(nrow=nrow(cont),ncol=9)
  for(i in 1:nrow(cont)){
    ###monolingual Russians
    mu[i,1] <-cont[i,]%*%rbind(1,0,0,0,0,0,inc[i],age)
    ###bilingual R/M Russians
    mu[i,2] <-cont[i,]%*%rbind(1,0,1,0,0,0,inc[i],age)
    ###bilingual R/U Russians
    mu[i,3] <-cont[i,]%*%rbind(1,0,0,1,0,0,inc[i],age)
    ###monolingual R Moldovans
    mu[i,4] <-cont[i,]%*%rbind(1,0,0,0,1,0,inc[i],age)
    ###bilingual R/M Moldovans
    mu[i,5] <-cont[i,]%*%rbind(1,0,1,0,1,0,inc[i],age)
    ###monolingual M Moldovans 
    mu[i,6] <-cont[i,]%*%rbind(1,1,1,0,1,0,inc[i],age)
    ###monolingual R Ukrainians
    mu[i,7] <-cont[i,]%*%rbind(1,0,0,0,1,1,inc[i],age)
    ###bilingual R/U Ukrainians
    mu[i,8] <-cont[i,]%*%rbind(1,0,0,1,1,1,inc[i],age)
    ###monolingual U Ukrainians
    mu[i,9] <-cont[i,]%*%rbind(1,1,0,1,1,1,inc[i],age)
    for(j in 1:9){
      pred[i,j]<-  1- pnorm(eps[i]  - mu[i,j])
    }
  }
  return(pred)
}

####Third threshold
eps<-as.matrix(cds[,166])

####Outcome coefficients
#####Supports integration with Moldova
contMol<- as.matrix(cds[,c(seq(44,103,10),seq(114,133,10))])
######Supports autonomy in Moldova
contAut<- as.matrix(cds[,c(seq(45,103,10),seq(115,133,10))])
####Supports confederation with Moldova
contCon<- as.matrix(cds[,c(seq(46,103,10),seq(116,133,10))])
########Supports independence
contInd<- as.matrix(cds[,c(seq(47,103,10),seq(117,133,10))])
##########Supports integration with Russia
contRus<- as.matrix(cds[,c(seq(48,103,10),seq(118,133,10))])
#######Supports autonomy within Russia
contRA<- as.matrix(cds[,c(seq(49,103,10),seq(119,133,10))])
########Supports confederation with Russia
contRC<- as.matrix(cds[,c(seq(50,103,10),seq(120,133,10))])
##########Supports integration with Ukraine
contUkr<- as.matrix(cds[,c(seq(51,103,10),seq(121,133,10))])
#######Supports autonomy within Ukraine
contUA<- as.matrix(cds[,c(seq(52,103,10),seq(122,133,10))])
########Supports confederation with Ukraine
contUC<- as.matrix(cds[,c(seq(53,103,10),seq(123,133,10))])

####Estimate for independence support; to estimate other outcomes use relevant matrix
mu<-bpreds(contInd,eps,inc,age)

####Estimate posterior median and HPD interval
mum<-apply(mu,2,median)
muhpd<-HPDinterval(as.mcmc(mu),prob=.9)

#########Dataset
Language<-as.factor(c("Russian monolingual","Russian/Moldovan bilingual","Russian/Ukrainian bilingual",
                      "Russian monolingual", "Russian/Moldovan bilingual", "Moldovan monolingual",
                      "Russian monolingual", "Russian/Ukrainian bilingual", "Ukrainian monolingual"))
Ethnicity <- as.factor(c(rep("Russian", 3),rep("Moldovan", 3),rep("Ukrainian", 3)))
mu<-data.frame(mum,muhpd,Language,Ethnicity)
names(mu)[1]<-"value"

mu$levels<-interaction(mu$Ethnicity,mu$Language)
mu$levels<-droplevels(mu$levels)
mu$levels = factor(mu$levels,levels(mu$levels)[c(3,7,6,4,8,9,2,5,1)])
levels(mu$levels)<-c("  Russian monolingual", "  Russian/Ukrainian bilingual", 
                     "  Russian/Moldovan bilingual", " Russian monolingual",        
                     " Russian/Ukrainian bilingual", " Ukrainian monolingual",      
                     "Russian monolingual", "Russian/Moldovan bilingual",  
                     "Moldovan monolingual")

mu$Ethnicity<-as.character(mu$Ethnicity)

####Figure M.1
ggplot(data = mu, aes(levels,value,color=Ethnicity, shape=Ethnicity))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(3.5,6.5)) + 
  theme_bw() + labs(x="", y="", color="",shape="") + 
  theme(text = element_text(size=15), legend.position="")  + 
  coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17)) 




###################################################
###############Latent variable analysis, Appendix N
###################################################

###Data file
hmod<-list(N=nrow(pmrdata), NDV=7, NDV2=2,t0=rep(0,4),T0=diag(.1,4), NL=3,
           y=array(c(as.matrix(pmrdata[,c("sr1","sr2","sr3","ur1","ur2","ur3","ur4")]),
                     as.matrix(pmrdata[,c("sm1","sm2","sm3","um1","um2","um3","um4")]),
                     as.matrix(pmrdata[,c("su1","su2","su3","uu1","uu2","uu3","uu4")])),
                   dim=c(nrow(pmrdata),7,3)),
           y2=array(c(as.matrix(pmrdata[,c("rsf","ruf")]),as.matrix(pmrdata[,c("msf","muf")]),
                      as.matrix(pmrdata[,c("usf","uuf")])),dim=c(nrow(pmrdata),2,3)))

####Initial threshold values
taustar=rbind(c(qnorm(mean(pmrdata$msf<2,na.rm=T)),qnorm(mean(pmrdata$msf<3,na.rm=T)),qnorm(mean(pmrdata$msf<4,na.rm=T)),qnorm(mean(pmrdata$msf<5,na.rm=T))),
              c(qnorm(mean(pmrdata$muf<2,na.rm=T)),qnorm(mean(pmrdata$muf<3,na.rm=T)),qnorm(mean(pmrdata$muf<4,na.rm=T)),qnorm(mean(pmrdata$muf<5,na.rm=T))))

##hmodel <- run.jags(method="parallel",model="lv_sp.txt", monitor=c("alpha","beta","tau","beta2","xi"), 
##                   inits=list(taustar=taustar), data=hmod, n.chains=4, burnin=10000, 
##                   sample=1250, thin=80, summarise=F,plots=F, modules=c("glm","lecuyer"))
##save(hmodel, file = "pf_sp.RData")
load("pf_sp.RData")

###Check convergence
cds<-as.mcmc.list(hmodel)
gd<-gelman.diag(cds[,1:24])
sum(gd[[1]][,1]>1.1)

####Convert to single chain
cds<-as.mcmc(hmodel)

##################Difficulty graphic
difficulty<-cds[,c(1:7,15:22)] 
difficulty[,8:15]<- -difficulty[,8:15]###Inverse order for thresholds

####Median and HPD interval
dm<-apply(difficulty,2,median)
dhpd<-HPDinterval(as.mcmc(difficulty))

###Data frame
d.ob<-data.frame(Question=c("Introductions","Talk about day","Discuss complex issues","Understand directions",
                            "Understand basic discussions","Understand quick speech","Understand hard TV", 
                            "Speaking T1", "Comprehension T1", "Speaking T2", "Comprehension T2", 
                            "Speaking T3", "Comprehension T3", "Speaking T4", "Comprehension T4"),dm,dhpd,Type=c(rep("Can do",7),rep(c("Speaking","Comprehension"),4)))
d.ob<-d.ob[order(d.ob[,2]),]
d.ob$Question<-factor(d.ob$Question,levels=d.ob[,1])

####Figure N.1.A
ggplot(data = d.ob, aes(Question, dm, color=Type))  + 
  geom_point(size = 2, position=position_dodge(width=0.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), size=1,show.legend = F,alpha=I(.5)) + 
  theme_bw() + labs(y="", x="") + scale_color_viridis_d() +
  theme(legend.position="top", text = element_text(size=20))  + coord_flip() + ggtitle("Posterior estimates of difficulty parameters")

##################Discrimination graphic
difficulty<-cds[,c(8:14,23:24)] 

####Median and HPDinterval
dm<-apply(difficulty,2,median)
dhpd<-HPDinterval(as.mcmc(difficulty))

###Data frame
d.ob<-data.frame(Question=c("Introductions","Talk about day","Discuss complex issues","Understand directions",
                            "Understand basic discussions","Understand quick speech","Understand hard TV","Speaking","Comprehension"),dm,dhpd)
d.ob<-d.ob[order(d.ob[,2]),]
d.ob$Question<-factor(d.ob$Question,levels=d.ob[,1])

####Figure N.1.A
ggplot(data = d.ob, aes(Question, dm))  + 
  geom_point(size = 2, position=position_dodge(width=0.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), size=1,show.legend = F,alpha=I(.25)) + 
  theme_bw() + labs(y="", x="") +
  theme(legend.position="top", text = element_text(size=20))  + coord_flip(ylim=c(0,14)) + ggtitle("Posterior estimates of discrimination \nparameters")


#######################################
###############Latent variable graphics
#######################################

######Russian estimates
xi<-cds[,25:601] 

####Median and HPD
xim<-apply(xi,2,median)
xihpd<-HPDinterval(as.mcmc(xi))

####Spoken proficiency
Speaking<-as.factor(pmrdata$rsf)

###Data frame
pl.ob<-data.frame(id=as.factor(seq(1:577)),xim,xihpd,Speaking)
pl.ob<-pl.ob[order(pl.ob[,2]),]
pl.ob$id<-factor(pl.ob$id,levels=pl.ob[,1])

####Figure N.2.A
ggplot(data = pl.ob, aes(id, xim,color=Speaking))  + 
  geom_point(size = 2, position=position_dodge(width=0.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), size=1,show.legend = F,alpha=I(.25)) + 
  theme_bw() + labs(y="", x="") + scale_fill_viridis_d() + scale_color_viridis_d() +
  theme(legend.position="top", text = element_text(size=20), axis.text.y = element_blank(),axis.ticks.y=element_blank())  + coord_flip(ylim=c(-2.5,2.5)) + ggtitle("Posterior estimates of Russian speaking proficiency")

###############Ukrainian estimates
xiU<-cds[,1179:1755] 

####Median and HPD
xiUm<-apply(xiU,2,median)
xiUhpd<-HPDinterval(as.mcmc(xiU))

####Spoken proficiency
Speaking<-as.factor(pmrdata$usf)

###Data frame
pl.ob<-data.frame(id=as.factor(seq(1:577)),xiUm,xiUhpd,Speaking)
pl.ob<-pl.ob[order(pl.ob[,2]),]
pl.ob$id<-factor(pl.ob$id,levels=pl.ob[,1])

####Figure N.2.A
ggplot(data = pl.ob, aes(id, xiUm,color=Speaking))  + 
  geom_point(size = 2, position=position_dodge(width=0.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), size=1,show.legend = F,alpha=I(.25)) + 
  theme_bw() + labs(y="", x="") + scale_fill_viridis_d() + scale_color_viridis_d() +
  theme(legend.position="top", text = element_text(size=20), axis.text.y = element_blank(),axis.ticks.y=element_blank())  + coord_flip(ylim=c(-2.5,2.5)) + ggtitle("Posterior estimates of Ukrainian speaking proficiency")

###############Moldovan estimates
xiM<-cds[,602:1178] 

####Median and HPD
xiMm<-apply(xiM,2,median)
xiMhpd<-HPDinterval(as.mcmc(xiM))

####Spoken proficiency
Speaking<-as.factor(pmrdata$msf)

###Data frame
pl.ob<-data.frame(id=as.factor(seq(1:577)),xiMm,xiMhpd,Speaking)
pl.ob<-pl.ob[order(pl.ob[,2]),]
pl.ob$id<-factor(pl.ob$id,levels=pl.ob[,1])

####Figure N.2.A
ggplot(data = pl.ob, aes(id, xiMm,color=Speaking))  + 
  geom_point(size = 2, position=position_dodge(width=0.5)) +
  geom_linerange(aes(ymin = lower, ymax = upper), size=1,show.legend = F,alpha=I(.25)) + 
  theme_bw() + labs(y="", x="") + scale_fill_viridis_d() + scale_color_viridis_d() +
  theme(legend.position="top", text = element_text(size=20), axis.text.y = element_blank(),axis.ticks.y=element_blank())  + coord_flip(ylim=c(-2.5,2.5)) + ggtitle("Posterior estimates of Moldovan speaking proficiency")

##############Add estimates to dataset
##pmrdata$xiR<-xim
##pmrdata$xiM<-xiMm
##pmrdata$xiU<-xiUm

##save(pmrdata,file="pmrdata.RData")

#########################################################################################
#############Regression analyses of support for separatism on latent linguistic abilities
#########################################################################################

#####Initial experimental threshold values
taustar<-polr(method="probit", as.factor(DV1) ~ xiR + xiM + xiU + NRussian + Ukrainian + Other + income + 
                age + male + highered + urban, data=pmrdata)$zeta

####Initial observational threshold values
gammastar<-matrix(nrow=10,ncol=4)
for(i in 1:10){
  gammastar[i,]<-polr(method="probit", as.factor(pmrdata[,i+23]) ~ pmrdata$xiR + pmrdata$xiM + 
                        pmrdata$xiU + pmrdata$NRussian + pmrdata$Ukrainian + pmrdata$Other + 
                        pmrdata$income + pmrdata$age + pmrdata$male + pmrdata$highered + 
                        pmrdata$urban)$zeta
}
gammastar<-apply(gammastar,2,mean)

###########Data file
hmod<-list(N=nrow(pmrdata), tg=pmrdata$tg, 
           y1=pmrdata$DV1, y2=as.matrix(pmrdata[,24:33]), 
           g0=rep(0,4), G0=diag(.01,4), 
           t0=rep(0,3), T0=diag(.01,3), 
           rf=pmrdata$xiR, mf=pmrdata$xiM, uf=pmrdata$xiU, NRussian=pmrdata$NRussian, 
           Other=pmrdata$Other, Ukrainian=pmrdata$Ukrainian, Male=pmrdata$male, 
           Age=pmrdata$age, Urban=pmrdata$urban,
           NIOTA=12, I0=diag(1,12), i0=rep(0,12), 
           NALPHA=7, A0=diag(1,7), a0=rep(0,7), 
           B0=diag(1,5), b0=rep(0,5), 
           HigherEd=pmrdata$highered, Married=pmrdata$married, Locality=pmrdata$locality, 
           Prof=pmrdata$prof, Resp=pmrdata$resp, House=pmrdata$house, ymin=pmrdata$ymin, 
           ymax=pmrdata$ymax, 
           bI0=rep(0,5), BI0=diag(1,5), ae1=rep(1,14))

##hmodel <- run.jags(method="parallel",model="pFc.txt", inits=list(gammastar=gammastar,
##                                                                 taustar=taustar), 
##                   monitor=c("alpha", "beta", "tau", "mu.a",
##                             "iota", "gamma",
##                             "pr", "pm",  "pu",
##                             "bI", "betaI","gammaI","alphaI",
##                             "pIe", "pIm","e1", "rIl","rIh",
##                             "Income"), data=hmod, n.chains=4, adapt=5000, burnin=50000, 
##                   sample=1000, thin=500, summarise=F,plots=F, modules=c("glm", "lecuyer"))
##save(hmodel, file = "pF1xi.RData")
load("pF1xi.RData")

########Analyze convergence
h1<-as.mcmc.list(hmodel)
gd<-gelman.diag(h1[,1:235],multivariate=F)
sum(gd[[1]][,1]>1.1)

####Convert to single chain
cds<-as.mcmc(hmodel)

###Average income and age across posterior draws
inc<- apply(as.matrix(cds[,233:809]),1,mean)
age<-mean(pmrdata$age)

############################experimental analysis
#############Function to estimate posterior probability respondent w different linguistic 
##and ethnic characteristics supports separatism across experimental conditions
#####Note that language is estimated at high and low modal values (.9, -1.15)
bpred<-function(cont,t1,t2,t3,eps,beta,inc,age){
  mu<-matrix(nrow=nrow(cont),ncol=36)
  pred<-matrix(nrow=nrow(cont),ncol=36)
  for(i in 1:nrow(cont)){
    ###monolingual Russians
    mu[i,1] <-cont[i,]%*%rbind(1,.9,-1.15,-1.15,0,0)
    mu[i,2] <-t1[i,]%*%rbind(1,.9,-1.15,-1.15,0,0)
    mu[i,3] <-t2[i,]%*%rbind(1,.9,-1.15,-1.15,0,0)
    mu[i,4] <-t3[i,]%*%rbind(1,.9,-1.15,-1.15,0,0)
    ###bilingual R/M Russians
    mu[i,5] <-cont[i,]%*%rbind(1,.9,.9,-1.15,0,0)
    mu[i,6] <-t1[i,]%*%rbind(1,.9,.9,-1.15,0,0)
    mu[i,7] <-t2[i,]%*%rbind(1,.9,.9,-1.15,0,0)
    mu[i,8] <-t3[i,]%*%rbind(1,.9,.9,-1.15,0,0)
    ###bilingual R/U Russians
    mu[i,9] <-cont[i,]%*%rbind(1,.9,-1.15,.9,0,0)
    mu[i,10] <-t1[i,]%*%rbind(1,.9,-1.15,.9,0,0)
    mu[i,11] <-t2[i,]%*%rbind(1,.9,-1.15,.9,0,0)
    mu[i,12] <-t3[i,]%*%rbind(1,.9,-1.15,.9,0,0)
    ###monolingual R Moldovans
    mu[i,13] <-cont[i,]%*%rbind(1,.9,-1.15,-1.15,1,0)
    mu[i,14] <-t1[i,]%*%rbind(1,.9,-1.15,-1.15,1,0)
    mu[i,15] <-t2[i,]%*%rbind(1,.9,-1.15,-1.15,1,0)
    mu[i,16] <-t3[i,]%*%rbind(1,.9,-1.15,-1.15,1,0)
    ###bilingual R/M Moldovans
    mu[i,17] <-cont[i,]%*%rbind(1,.9,.9,-1.15,1,0)
    mu[i,18] <-t1[i,]%*%rbind(1,.9,.9,-1.15,1,0)
    mu[i,19] <-t2[i,]%*%rbind(1,.9,.9,-1.15,1,0)
    mu[i,20] <-t3[i,]%*%rbind(1,.9,.9,-1.15,1,0)
    ###monolingual M Moldovans 
    mu[i,21] <-cont[i,]%*%rbind(1,-1.15,.9,-1.15,1,0)
    mu[i,22] <-t1[i,]%*%rbind(1,-1.15,.9,-1.15,1,0)
    mu[i,23] <-t2[i,]%*%rbind(1,-1.15,.9,-1.15,1,0)
    mu[i,24] <-t3[i,]%*%rbind(1,-1.15,.9,-1.15,1,0)
    ###monolingual R Ukrainians
    mu[i,25] <-cont[i,]%*%rbind(1,.9,-1.15,-1.15,1,1)
    mu[i,26] <-t1[i,]%*%rbind(1,.9,-1.15,-1.15,1,1)
    mu[i,27] <-t2[i,]%*%rbind(1,.9,-1.15,-1.15,1,1)
    mu[i,28] <-t3[i,]%*%rbind(1,.9,-1.15,-1.15,1,1)
    ###bilingual R/U Ukrainians
    mu[i,29] <-cont[i,]%*%rbind(1,.9,-1.15,.9,1,1)
    mu[i,30] <-t1[i,]%*%rbind(1,.9,-1.15,.9,1,1)
    mu[i,31] <-t2[i,]%*%rbind(1,.9,-1.15,.9,1,1)
    mu[i,32] <-t3[i,]%*%rbind(1,.9,-1.15,.9,1,1)
    ###monolingual U Ukrainians
    mu[i,33] <-cont[i,]%*%rbind(1,-1.15,-1.15,.9,1,1)
    mu[i,34] <-t1[i,]%*%rbind(1,-1.15,-1.15,.9,1,1)
    mu[i,35] <-t2[i,]%*%rbind(1,-1.15,-1.15,.9,1,1)
    mu[i,36] <-t3[i,]%*%rbind(1,-1.15,-1.15,.9,1,1)
    for(j in 1:36){
      pred[i,j]<-  1- pnorm(eps[i]  - (mu[i,j] + beta[i,]%*%rbind(inc[i],age)))
    }
  }
  return(pred)
}

####Third threshold
eps<-as.matrix(cds[,36])

#####Experimental coefficients
cont<- as.matrix(cds[,seq(1,24,4)])
t1<- as.matrix(cds[,seq(2,24,4)])
t2<- as.matrix(cds[,seq(3,24,4)])
t3<- as.matrix(cds[,seq(4,24,4)])
###Treatment invariant coefficients
beta<- as.matrix(cds[,29:30])

#####Estimate posterior probability of support for separatism
mu<-bpred(cont,t1,t2,t3,eps,beta,inc,age)

#####Estimate median and HPD interval
mum<-apply(mu,2,median)
muhpd<-HPDinterval(as.mcmc(mu),prob=.95)

#####Data frame
Condition<-rep(c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian"),9)
Ethnicity<-c(rep("Russian",12),rep("Moldovan",12),rep("Ukrainian",12))
Language<-factor(c(rep("Russian monolingual",4),rep("Russian/Moldovan bilingual",4),
            rep("Russian/Ukrainian bilingual",4), rep("Russian monolingual",4),
            rep("Russian/Moldovan bilingual",4), rep("Moldovan monolingual",4),
            rep("Russian monolingual",4),rep("Russian/Ukrainian bilingual",4),
            rep("Ukrainian monolingual",4)))
mu<-data.frame(mum,muhpd,Condition, Ethnicity,Language)
names(mu)[1]<-"value"

####Control
muc<-mu[mu$Condition=="Control",]
muc$levels<-interaction(muc$Ethnicity,muc$Language)
muc$levels<-droplevels(muc$levels)
muc$levels = factor(muc$levels,levels(muc$levels)[c(3,7,6,4,8,9,2,5,1)])
levels(muc$levels)<-c("  Russian monolingual", "  Russian/Ukrainian bilingual", 
                      "  Russian/Moldovan bilingual", " Russian monolingual",        
                      " Russian/Ukrainian bilingual", " Ukrainian monolingual",      
                      "Russian monolingual", "Russian/Moldovan bilingual",  
                      "Moldovan monolingual")

muc$Ethnicity<-as.character(muc$Ethnicity)

###N.4.A
ggplot(data = muc, aes(levels,value,color=Ethnicity, shape=Ethnicity))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(3.5,6.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="") + 
  theme(text = element_text(size=15), legend.position=c(.25,.25))  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

###Subgroup
mu$levels<-interaction(mu$Ethnicity,mu$Language,mu$Condition)
mu$levels<-droplevels(mu$levels)
mu$levels = factor(mu$levels,sort(levels(mu$levels)))
mu$levels = factor(mu$levels,levels(mu$levels)[c(13:16,21:24,17:20,25:36,5:12,1:4)])

levels(mu$levels)<-c("Control", "T_Moldovan", "T_Russian", "T_Ukrainian", 
                     " Control", " T_Moldovan", " T_Russian", " T_Ukrainian",
                     "  Control", "  T_Moldovan", "  T_Russian", "  T_Ukrainian",
                     "   Control", "   T_Moldovan", "   T_Russian", "   T_Ukrainian",
                     "    Control", "    T_Moldovan", "    T_Russian", "    T_Ukrainian",
                     "     Control", "     T_Moldovan", "     T_Russian", "     T_Ukrainian",
                     "      Control", "      T_Moldovan", "      T_Russian", "      T_Ukrainian",
                     "       Control", "       T_Moldovan", "       T_Russian", "       T_Ukrainian",
                     "        Control", "        T_Moldovan", "        T_Russian", "        T_Ukrainian")

mu$Ethnicity<-as.character(mu$Ethnicity)

####Figure N.4.B
ggplot(data = mu[which(mu$Ethnicity=="Russian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))


###Figure N.4.C
ggplot(data = mu[which(mu$Ethnicity=="Moldovan"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))

###Figure N.4.D
ggplot(data = mu[which(mu$Ethnicity=="Ukrainian"),], 
       aes(levels,value,shape=Language, color=Language))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(4.5,8.5)) + 
  theme_bw() + labs(y="", x="", color="",shape="",linetype="") + guides(colour = guide_legend(nrow = 2)) + 
  theme(text = element_text(size=15), legend.position="top")  + coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17))


#######################Observational analyses
#Function to estimate posterior probability respondent w different ethnic and linguistic
#characteristics supports a given separatist outcome
bpreds<-function(cont,eps,inc,age){
  mu<-matrix(nrow=nrow(cont),ncol=9)
  pred<-matrix(nrow=nrow(cont),ncol=9)
  for(i in 1:nrow(cont)){
    ###monolingual Russians
    mu[i,1] <-cont[i,]%*%rbind(1,.9,-1.15,-1.15,0,0,inc[i],age)
    ###bilingual R/M Russians
    mu[i,2] <-cont[i,]%*%rbind(1,.9,.9,-1.15,0,0,inc[i],age)
    ###bilingual R/U Russians
    mu[i,3] <-cont[i,]%*%rbind(1,.9,-1.15,.9,0,0,inc[i],age)
    ###monolingual R Moldovans
    mu[i,4] <-cont[i,]%*%rbind(1,.9,-1.15,-1.15,1,0,inc[i],age)
    ###bilingual R/M Moldovans
    mu[i,5] <-cont[i,]%*%rbind(1,.9,.9,-1.15,1,0,inc[i],age)
    ###monolingual M Moldovans 
    mu[i,6] <-cont[i,]%*%rbind(1,-1.15,.9,-1.15,1,0,inc[i],age)
    ###monolingual R Ukrainians
    mu[i,7] <-cont[i,]%*%rbind(1,.9,-1.15,-1.15,1,1,inc[i],age)
    ###bilingual R/U Ukrainians
    mu[i,8] <-cont[i,]%*%rbind(1,.9,-1.15,.9,1,1,inc[i],age)
    ###monolingual U Ukrainians
    mu[i,9] <-cont[i,]%*%rbind(1,-1.15,-1.15,.9,1,1,inc[i],age)
    for(j in 1:9){
      pred[i,j]<-  1- pnorm(eps[i]  - mu[i,j])
    }
  }
  return(pred)
}

####Third threshold
eps<-as.matrix(cds[,166])
###Support integration w Moldova
contMol<- as.matrix(cds[,c(seq(44,103,10),seq(114,133,10))])
###Support autonomy in Moldova
contAut<- as.matrix(cds[,c(seq(45,103,10),seq(115,133,10))])
###Support confederation w Moldova
contCon<- as.matrix(cds[,c(seq(46,103,10),seq(116,133,10))])
###Support independence
contInd<- as.matrix(cds[,c(seq(47,103,10),seq(117,133,10))])
###Support integration w Russia
contRus<- as.matrix(cds[,c(seq(48,103,10),seq(118,133,10))])
###Support autonomy in Russia
contRA<- as.matrix(cds[,c(seq(49,103,10),seq(119,133,10))])
###Support confederation w Russia
contRC<- as.matrix(cds[,c(seq(50,103,10),seq(120,133,10))])
###Support integration w Ukraine
contUkr<- as.matrix(cds[,c(seq(51,103,10),seq(121,133,10))])
###Support autonomy in Ukraine
contUA<- as.matrix(cds[,c(seq(52,103,10),seq(122,133,10))])
###Support confederation w Ukraine
contUC<- as.matrix(cds[,c(seq(53,103,10),seq(123,133,10))])


######Estimate probability of support for independence. Replace with relevant matrix to estimate
#other outcomes
mu<-bpreds(contInd,eps,inc,age)

####Estimate median and HPD interval
mum<-apply(mu,2,median)
muhpd<-HPDinterval(as.mcmc(mu),prob=.95)

####Data frame
Language<-as.factor(c("Russian monolingual","Russian/Moldovan bilingual","Russian/Ukrainian bilingual",
                      "Russian monolingual", "Russian/Moldovan bilingual", "Moldovan monolingual",
                      "Russian monolingual", "Russian/Ukrainian bilingual", "Ukrainian monolingual"))

Ethnicity <- as.factor(c(rep("Russian", 3),rep("Moldovan", 3),rep("Ukrainian", 3)))
mu<-data.frame(mum,muhpd,Language,Ethnicity)
names(mu)[1]<-"value"

mu$levels<-interaction(mu$Ethnicity,mu$Language)
mu$levels<-droplevels(mu$levels)
mu$levels = factor(mu$levels,levels(mu$levels)[c(3,7,6,4,8,9,2,5,1)])
levels(mu$levels)<-c("  Russian monolingual", "  Russian/Ukrainian bilingual", 
                     "  Russian/Moldovan bilingual", " Russian monolingual",        
                     " Russian/Ukrainian bilingual", " Ukrainian monolingual",      
                     "Russian monolingual", "Russian/Moldovan bilingual",  
                     "Moldovan monolingual")
mu$Ethnicity<-as.character(mu$Ethnicity)

####Figure N.3
ggplot(data = mu, aes(levels,value,color=Ethnicity, shape=Ethnicity))  + 
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), size=1,width = .01, 
                position=position_dodge(width=0.5)) + 
  geom_vline(xintercept = c(3.5,6.5)) + 
  theme_bw() + labs(x="", y="", color="",shape="") + 
  theme(text = element_text(size=15), legend.position="")  + 
  coord_flip(ylim=c(0,1)) + scale_color_grey() +
  scale_shape_manual(values=c(15,16,17)) 